Tissue fusion and enhanced genotypic diversity supports survival of juvenile Pocillopora acuta corals under thermal stress

Experimental Design

Survival of juvenile corals is vital for the recovery and resilience of reefs as global climate change and local stressors threaten the persistence of these ecosystems. However, newly settled coral recruits are at high risk for mortality due to their small size. Tissue fusion with neighboring recruits provides a pathway to rapidly increasing total colony size faster than would be achieved by growth alone. As self-recognition systems do not fully develop until several months post-settlement, fusion between newly settled juvenile colonies provides an opportunity for an increase in genetic diversity along with an increase in size. Although fusion provides advantages for survival across several species, it is unknown whether an increase in size or genetic diversity drives increases in survival. As ocean temperatures continue to increase, these potential benefits of fusion may provide ecological advantages during thermal stress events. In this study, We test these concepts by investigating the implications of fusion between newly settled Pocillopora acuta juveniles. We manipulated the genotypic richness and number of juveniles involved in tissue fusion to analyze influences on survival and growth in thermal stress. In ambient conditions, tissue fusion provided a significant survival advantage regardless of genotypic diversity, while negative competitive effects existed between closely settled, but non-fused, juveniles. In the presence of thermal stress, tissue fusion significantly delayed juvenile mortality and fusion between multiple genotypes provided an additional survival advantage. Delayed mortality as a result by tissue fusion and elevated genotypic diversity provides a window of opportunity for increased survival of fused juveniles. These results indicate that tissue fusion provides important ecological advantages in thermal stress for P. acuta, which are further enhanced by increased genetic diversity. Further research is required to determine energetic and genetic mechanisms that contribute to increases in survivorship as a result of tissue fusion.

Please see manuscript for detailed methodology

In this experiment, settled juvenile corals were secured on slides either as solitary individuals or in groups. Groups were composed of 2, 3, or 4 juveniles either from the same genotype or of multiple genotypes. See genetic analyzes for distinction of genotypes. These slides were then placed in an indoor tank system with light and flow-through seawater at ambient temperature for a period of 15 days, allowing them to grow and have the opportunity to fuse. This period is referred to as the “grow out period”. At the end of the grow-out period, we measured survivorship (alive or dead), growth (as number of polyps per juvenile), and fusion rate (fusion success or no fusion). Fusion was identified by continuous pigmentation across boundaries between juveniles with clear tissue connectivitiy. We then analyze the effect of genotypic richness and fusion on survival and growth. All responses in this experiment were measured at the level of the individual juvenile, not at the level of the slide community.

After the 15-day grow-out period, we exposed juveniles on slides to either a continuation of the ambient treatment (daily maximum temperature = 27.4 °C) or exposure to a high temperature (daily maximum temperature = 30.8 °C). At multiple points during this “stress period”, we measured survivorship (alive or dead), growth (numer of polyps per juvenile), and fusion persistence (fusion success or no fusion). We analyzed the effect of temperature, fusion (classified as the fusion status at the start of the stress period), and genotypic richness on survivorship and growth. Additionally, we analyzed the effect of temperature and genotypic richness on fusion persistence.

Experimental design schematic. Circles represent individual juveniles and colors represent genotypes.

Example of juveniles settled on slides as a group.

Slides in an indoor tank during the grow-out period.

Data Analysis

Analysis Setup

Set up workspace, set options, and load required packages.

rm(list=ls(all=TRUE)) 

library(gridExtra)
library(multcomp)
library(emmeans)
library(ggplot2)
library(dplyr)
library(effects)
library(plyr)
library(lattice)
library(glmmTMB)
library(effects)
library(lsmeans)
library(MuMIn)
library(MASS)
library(car)
library(FSA)
library(lmerTest)
library(lme4)
library(blmeco)
library(ggsci)
library(coxme)
library(survival)
library(cowplot)
library(RColorBrewer)
library(tidyverse)
library(partitionBEFsp)
library(factoextra)
library(ggfortify)
library(cowplot)
library(survminer)

1. Grow-Out Period: Survivorship

Load data

Survivorship is measured as binary (alive = 1, dead = 0) on each individual spat in the experiment.

#load data
rearing<-read.csv("Data/GrowOut_Responses_PostGenotyping.csv", header=TRUE, sep=",", na.strings="NA") #load data
rearing$Quadrant<-as.factor(rearing$Quadrant) #change to factors
rearing$Diversity.Type<-relevel(rearing$Diversity.Type, ref="Individual")
rearing$Community<-relevel(rearing$Community, ref="Individual")

Subset the final timepoint for additional analyses.

final<-rearing[which(rearing$Timepoint == "Final"),] #subset final timepoint
group<-rearing[which(rearing$Community == "Group"),] #subset juveniles in "groups"
group_final<-group[which(group$Timepoint == "Final"),] #subset the final timepoint for groups

Analysis

Analyze survival at the end (final timepoint) of the grow-out period with a logistic binomial regression model.

Main effects: Genotypic richness, community size (either individual juveniles or groups of juveniles), fusion partners (the number of juveniles each juvenile was fused with), and the interaction of genotypic richness and fusion.

Random effects: Parent colony nested within parent colony collection reef, tank, and slide.

Here, generate the summary of the model and the anova table output for use of p-values.

#with final dataset
survmod2<-glmer(Survivorship ~ Richness + Community + Fusion.Partners + Richness:Fusion.Partners + (1|Parent.Site/Colony) + (1|Tank) + (1|Slide), data=final, family=binomial)
summary(survmod2)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: 
## Survivorship ~ Richness + Community + Fusion.Partners + Richness:Fusion.Partners +  
##     (1 | Parent.Site/Colony) + (1 | Tank) + (1 | Slide)
##    Data: final
## 
##      AIC      BIC   logLik deviance df.resid 
##    224.6    255.2   -103.3    206.6      212 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -10.6985  -0.6214   0.3067   0.5208   1.4376 
## 
## Random effects:
##  Groups             Name        Variance  Std.Dev.
##  Slide              (Intercept) 3.391e-01 0.582352
##  Colony:Parent.Site (Intercept) 4.784e-03 0.069167
##  Parent.Site        (Intercept) 3.560e-06 0.001887
##  Tank               (Intercept) 1.537e-01 0.392056
## Number of obs: 221, groups:  
## Slide, 115; Colony:Parent.Site, 21; Parent.Site, 6; Tank, 5
## 
## Fixed effects:
##                          Estimate Std. Error z value Pr(>|z|)   
## (Intercept)                1.3942     0.5005   2.786  0.00534 **
## Richness                   0.0457     0.2590   0.176  0.85997   
## CommunityGroup            -1.9190     0.6414  -2.992  0.00277 **
## Fusion.Partners            3.3668     1.0910   3.086  0.00203 **
## Richness:Fusion.Partners  -0.5213     0.3525  -1.479  0.13924   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Rchnss CmmntG Fsn.Pr
## Richness    -0.477                     
## CommuntyGrp -0.207 -0.608              
## Fusn.Prtnrs -0.114  0.411 -0.482       
## Rchnss:Fs.P  0.153 -0.458  0.393 -0.929
## convergence code: 0
## Model failed to converge with max|grad| = 0.0112085 (tol = 0.001, component 1)

Anova(survmod2)

## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: Survivorship
##                            Chisq Df Pr(>Chisq)    
## Richness                  0.3173  1   0.573223    
## Community                 8.9516  1   0.002772 ** 
## Fusion.Partners          21.3756  1  3.775e-06 ***
## Richness:Fusion.Partners  2.1863  1   0.139245    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

There is a significant difference between individual and groups of juveniles with a significant effect of fusion on survival.

Display effect estimates for 1) Richness and Fusion Partners; 2) Fusion Partners; 3) Individual vs Groups

#final dataset
summary(effect(c("Richness", "Fusion.Partners"), xlevels=4, survmod2))

## 
##  Richness*Fusion.Partners effect
##         Fusion.Partners
## Richness         0         1         2         3
##        1 0.5061314 0.9463480 0.9967168 0.9998087
##        2 0.5175489 0.9164084 0.9911532 0.9991274
##        3 0.5289481 0.8720208 0.9763849 0.9960299
##        4 0.5403171 0.8089765 0.9384931 0.9821349
## 
##  Lower 95 Percent Confidence Limits
##         Fusion.Partners
## Richness         0         1         2         3
##        1 0.3130708 0.8186676 0.9501058 0.9867223
##        2 0.3661859 0.8084490 0.9462577 0.9856340
##        3 0.3374863 0.7366122 0.9027604 0.9645174
##        4 0.2661940 0.5541801 0.6689295 0.7279225
## 
##  Upper 95 Percent Confidence Limits
##         Fusion.Partners
## Richness         0         1         2         3
##        1 0.6973806 0.9856964 0.9997934 0.9999973
##        2 0.6657585 0.9660745 0.9985992 0.9999477
##        3 0.7122548 0.9431848 0.9945985 0.9995683
##        4 0.7920387 0.9351828 0.9913962 0.9991155

summary(effect("Fusion.Partners", xlevels=4, survmod2))

## 
##  Fusion.Partners effect
## Fusion.Partners
##         0         1         2         3 
## 0.5169809 0.9182039 0.9915778 0.9991908 
## 
##  Lower 95 Percent Confidence Limits
## Fusion.Partners
##         0         1         2         3 
## 0.3653903 0.8098959 0.9468856 0.9858654 
## 
##  Upper 95 Percent Confidence Limits
## Fusion.Partners
##         0         1         2         3 
## 0.6655093 0.9672972 0.9987155 0.9999543

summary(effect("Community", xlevels=4, survmod2))

## 
##  Community effect
## Community
## Individual      Group 
##  0.9532330  0.7494544 
## 
##  Lower 95 Percent Confidence Limits
## Community
## Individual      Group 
##  0.8445862  0.5984934 
## 
##  Upper 95 Percent Confidence Limits
## Community
## Individual      Group 
##  0.9870881  0.8571990

Conduct diagnostics. Passes test for overdispersion (value = 0.9).

dispersion_glmer(survmod2) #no evidence of overdispersion
plot(residuals(survmod2)) + abline(h=0, lty=2) #Dispersed randomly

Plotting

Create plot for effect of community type (individual juveniles vs. groups of juveniles).

eff.surv1 <- predictorEffect("Community", survmod2) #set x axis groups
surv.effects1<-plot(eff.surv1,
                   lines=list(multiline=TRUE, #color and type of lines
                              col=c("black"), 
                              lty=1), 
                   confint=list(style="bars", width=4, col="black"), #confidence interval style
                   lwd=4,
                   axes=list(y=list(lim=c(0, 1))),
                   type="response", #scale of response - "response" gives in original data scale
                   ylab=expression(bold("Prob(Survivorship)")), 
                   legend.position="right",
                   xlab=expression(bold("Number of Juveniles")), 
                   main=""); surv.effects1

Create plot to provide direct comparison of survival between individual juveniles, fused juveniles, and unfused juveniles. Display summary table of means.

#add categories to identify individuals, fused, or nonfused 
final$Fusion.Type<-ifelse(final$Fusion %in% 1, "Fused",
                   ifelse(final$Fusion %in% 0, "Unfused", NA))

final$Fusion.Type<-ifelse(final$Community %in% c("Individual"), "Individual", final$Fusion.Type)
  
#generate summary table 
grow_plot_surv <- ddply(final, c("Fusion.Type"), summarise,
                   N    = length(Survivorship[!is.na(Survivorship)]),
                   mean = mean(Survivorship, na.rm=TRUE),
                   sd   = sd(Survivorship, na.rm=TRUE),
                   se   = sd / sqrt(N)) 
                
grow_plot_surv

##   Fusion.Type  N      mean        sd         se
## 1       Fused 88 0.9431818 0.2328215 0.02481886
## 2  Individual 58 0.7758621 0.4206554 0.05523477
## 3     Unfused 75 0.3466667 0.4791133 0.05532324

#plot survival 
mean_growout_plot<-ggplot(data=grow_plot_surv, aes(x=Fusion.Type, y=mean)) +
  geom_point(size=6, position=position_dodge(0.03), color="black") +
  geom_errorbar(aes(ymin=mean-se, ymax=mean+se), width=0, linetype="solid", position=position_dodge(0.03), size=2) +
  theme_classic()+
  theme(text = element_text(size = 18, color="black"))+
  theme(axis.text = element_text(size = 18, color="black"))+
  theme(legend.position = "none")+
  theme(axis.title = element_text(size = 22, color="black", face="bold"))+
  ylab(expression(bold(paste("Survivorship")))) +
  theme(legend.title=element_blank())+
  theme(legend.text = element_text(size=18))+
  theme(plot.margin=unit(c(1,0.1,0,0.1), "cm"))+
  scale_y_continuous(limits=c(0,1),breaks=c(0.0,0.2,0.4,0.6,0.8,1.0))+
  scale_x_discrete(limits=c("Individual", "Fused", "Unfused"))+
  theme(axis.title.y = element_text(margin = margin(t = 0, r = 7, b = 0, l = 1))) +
  theme(axis.title.x = element_text(margin = margin(t = 10, r = 0, b = 0, l = 0))) +
  xlab("Juvenile Type"); mean_growout_plot

ggsave("Figures/grow_out_means.pdf", plot=mean_growout_plot, height=11, width=6, units = c("in"), dpi=500) #output figure

Create plot for the effect of fusion partners on survival within groups of juveniles.

survmod2b<-glmer(Survivorship ~ Fusion.Partners + (1|Parent.Site/Colony) + (1|Tank) + (1|Slide), data=final, family=binomial, subset=c(Community=="Group"))

eff.surv2b <- predictorEffect("Fusion.Partners", survmod2b, xlevels=4, rug=TRUE)
surv.effects2b<-plot(eff.surv2b,
                   lines=list(multiline=TRUE, 
                              col=c("black"), 
                              lty=1), 
                   confint=list(style="bands", alpha=0.1), 
                   lwd=4,
                   axes=list(y=list(lim=c(0, 1))),
                   type="response", 
                   ylab=expression(bold("")), 
                   legend.position="right",
                   xlab=expression(bold("Fusion Partners")), 
                   main="",
                   lattice=list(key.args=list(space="none",
                                              border=FALSE, 
                                              title=expression(bold("")),
                                              cex=1, 
                                              cex.title=1)));surv.effects2b

Within groups, there was lowest survival in juveniles without fusion partners (unfused juveniles). There was a survival advantage of juveniles fused with others.

Generate mean values of survivorship for 1) within “group” juveniles, mean values for number of fusion partners; and 2) within “group” juveniles, mean values for genotypic richnes.

meansurv <- ddply(final, c("Community"), summarise,
                   N    = length(Survivorship[!is.na(Survivorship)]),
                   mean = mean(Survivorship, na.rm=TRUE),
                   sd   = sd(Survivorship, na.rm=TRUE),
                   se   = sd / sqrt(N), 
                   max = max(Survivorship, na.rm=TRUE),
                   lower = mean-sd, 
                   upper = mean+sd
                   )

meansurv2 <- ddply(group_final, c("Fusion.Partners"), summarise,
                   N    = length(Survivorship[!is.na(Survivorship)]),
                   mean = mean(Survivorship, na.rm=TRUE),
                   sd   = sd(Survivorship, na.rm=TRUE),
                   se   = sd / sqrt(N), 
                   max = max(Survivorship, na.rm=TRUE),
                   lower = mean-sd, 
                   upper = mean+sd
                   )
meansurv2

##   Fusion.Partners  N      mean        sd         se max      lower    upper
## 1               0 75 0.3466667 0.4791133 0.05532324   1 -0.1324466 0.825780
## 2               1 48 0.9583333 0.2019409 0.02914766   1  0.7563924 1.160274
## 3               2 24 0.9166667 0.2823299 0.05763034   1  0.6343368 1.198997
## 4               3 16 0.9375000 0.2500000 0.06250000   1  0.6875000 1.187500

meansurv3 <- ddply(group_final, c("Richness"), summarise,
                   N    = length(Survivorship[!is.na(Survivorship)]),
                   mean = mean(Survivorship, na.rm=TRUE),
                   sd   = sd(Survivorship, na.rm=TRUE),
                   se   = sd / sqrt(N), 
                   max = max(Survivorship, na.rm=TRUE),
                   lower = mean-sd, 
                   upper = mean+sd
                   )
meansurv3

##   Richness  N      mean        sd         se max      lower    upper
## 1        1 57 0.6842105 0.4689614 0.06211546   1 0.21524911 1.153172
## 2        2 26 0.6153846 0.4961389 0.09730085   1 0.11924568 1.111524
## 3        3 56 0.7142857 0.4558423 0.06091449   1 0.25844341 1.170128
## 4        4 24 0.5833333 0.5036102 0.10279899   1 0.07972318 1.086943

2. Grow-Out Period: Fusion Rate

Fusion rate is analyzed as a binary response with 1 indicating fused and 0 indicating not fused.

Analysis

Analyze fusion rate at the end of the grow-out period with a logistic binomial regression model. Here, we analyze fusion as a product of genotypic richness within “group” juveniles.

Main effects: Genoyptic richness
Random effects: Parent colony nested within parent collection site, tank, slide

Display model summary and anova table for use of p-values.

#with final dataset
fusemod1<-glmer(Fusion ~ Richness + (1|Parent.Site/Colony) + (1|Tank) + (1|Slide), data=final, family=binomial, subset=c(Community=="Group"))
summary(fusemod1)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: Fusion ~ Richness + (1 | Parent.Site/Colony) + (1 | Tank) + (1 |  
##     Slide)
##    Data: final
##  Subset: c(Community == "Group")
## 
##      AIC      BIC   logLik deviance df.resid 
##    184.8    203.4    -86.4    172.8      157 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.6021 -0.2575  0.1468  0.2860  1.1850 
## 
## Random effects:
##  Groups             Name        Variance Std.Dev.
##  Slide              (Intercept) 14.8391  3.8522  
##  Colony:Parent.Site (Intercept)  0.4334  0.6583  
##  Parent.Site        (Intercept)  1.0245  1.0122  
##  Tank               (Intercept)  2.7852  1.6689  
## Number of obs: 163, groups:  
## Slide, 57; Colony:Parent.Site, 21; Parent.Site, 6; Tank, 5
## 
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)
## (Intercept)   1.7587     2.1221   0.829    0.407
## Richness     -0.3583     0.7338  -0.488    0.625
## 
## Correlation of Fixed Effects:
##          (Intr)
## Richness -0.832

Anova(fusemod1, type="II")

## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: Fusion
##           Chisq Df Pr(>Chisq)
## Richness 0.2384  1     0.6253

There is no effect of genotypic richness on fusion success.

Conduct diagnostics. There is no evidence of overdispersion.

dispersion_glmer(fusemod1) #no overdispersion, less than 1.4
plot(residuals(fusemod1)) + abline(h=0, lty=2)

Plotting

There is no effect of genotypic diversity on fusion rates. Plot this relationship.

eff.fuse1 <- predictorEffect("Richness", fusemod1, xlevels=4, rug=TRUE) #set x axis
fuse.effects1<-plot(eff.fuse1,
                   lines=list(multiline=TRUE, #color lines
                              col=c("black"), 
                              lty=1),
                   confint=list(style="bands", alpha=0.1), #set conf int style
                   lwd=4,
                   type="response", #set resposne scale
                   ylab=expression(bold("Prob(Successful Fusion)")), 
                   xlab=expression(bold("Richness")), 
                   main="",
                   lattice=list(key.args=list(space="right",
                                              border=FALSE, 
                                              title=expression(bold("Genotypic \nRichness")),
                                              cex=1, 
                                              cex.title=1)));fuse.effects1

Generate a summary of proportion fused at the end of the grow-out period by genotypic richness.

meanfusion <- ddply(group_final, c("Richness"), summarise,
                   N    = length(Fusion[!is.na(Fusion)]),
                   mean = mean(Fusion, na.rm=TRUE),
                   sd   = sd(Fusion, na.rm=TRUE),
                   se   = sd / sqrt(N), 
                   max = max(Fusion, na.rm=TRUE),
                   lower = mean-sd, 
                   upper = mean+sd
                   )
meanfusion

##   Richness  N      mean        sd         se max       lower    upper
## 1        1 57 0.5263158 0.5037454 0.06672270   1  0.02257042 1.030061
## 2        2 26 0.6153846 0.4961389 0.09730085   1  0.11924568 1.111524
## 3        3 56 0.5357143 0.5032363 0.06724778   1  0.03247801 1.038951
## 4        4 24 0.5000000 0.5107539 0.10425721   1 -0.01075392 1.010754

Fusion occured in approximately 50% of juveniles with no effect of genotypic richness.

3. Grow-Out Period: Polyp Expansion (Growth)

Analysis

Growth is measured as polyp expansion in each spat by the number of polyps in each spat by the end of the grow-out period. As all juveniels started as one primary polyp, there is no need to correct for initial size.

Analyze polyp expansion with a poisson mixed effects model in glmer.

Main effects: Genotypic richness, community size (individuals vs groups), and number of fusion partners.
Random effects: Parent colony nested with parent collection site, tank, and slide.

Display model summary and anova table for use of p-values.

polypmod1<-glmer(Polyps ~ Richness + Community + Fusion.Partners + Richness:Fusion.Partners + (1|Parent.Site/Colony) + (1|Tank) + (1|Slide), data=final, family = poisson)
summary(polypmod1)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: 
## Polyps ~ Richness + Community + Fusion.Partners + Richness:Fusion.Partners +  
##     (1 | Parent.Site/Colony) + (1 | Tank) + (1 | Slide)
##    Data: final
## 
##      AIC      BIC   logLik deviance df.resid 
##    590.7    618.2   -286.4    572.7      148 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.10224 -0.55740 -0.03806  0.43006  1.42031 
## 
## Random effects:
##  Groups             Name        Variance  Std.Dev. 
##  Slide              (Intercept) 4.527e-09 6.728e-05
##  Colony:Parent.Site (Intercept) 7.002e-02 2.646e-01
##  Parent.Site        (Intercept) 5.196e-03 7.208e-02
##  Tank               (Intercept) 7.883e-03 8.879e-02
## Number of obs: 157, groups:  
## Slide, 92; Colony:Parent.Site, 21; Parent.Site, 6; Tank, 5
## 
## Fixed effects:
##                           Estimate Std. Error z value Pr(>|z|)    
## (Intercept)               1.245663   0.150290   8.288   <2e-16 ***
## Richness                  0.002656   0.089657   0.030    0.976    
## CommunityGroup           -0.089969   0.174884  -0.514    0.607    
## Fusion.Partners          -0.013179   0.168666  -0.078    0.938    
## Richness:Fusion.Partners  0.014306   0.054978   0.260    0.795    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Rchnss CmmntG Fsn.Pr
## Richness    -0.592                     
## CommuntyGrp  0.120 -0.689              
## Fusn.Prtnrs -0.353  0.642 -0.711       
## Rchnss:Fs.P  0.429 -0.768  0.656 -0.936
## convergence code: 0
## boundary (singular) fit: see ?isSingular

Anova(polypmod1, type="II")

## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: Polyps
##                           Chisq Df Pr(>Chisq)
## Richness                 0.1281  1     0.7204
## Community                0.2647  1     0.6069
## Fusion.Partners          0.2204  1     0.6387
## Richness:Fusion.Partners 0.0677  1     0.7947

There is no effect of genotypic richness and fusion and no difference in growth between individual or groups of juveniles.

Conduct diagnostics with a qq-plot to assess normality. Normality is met.

qqPlot(residuals(polypmod1))

## 161 165 
##  32  33

Display mean values of polyp expansion at the end of the grow-out period.

meangrowth <- ddply(final, c("Community"), summarise,
                   N    = length(Polyps[!is.na(Polyps)]),
                   mean = mean(Polyps, na.rm=TRUE),
                   sd   = sd(Polyps, na.rm=TRUE),
                   se   = sd / sqrt(N), 
                   max = max(Polyps, na.rm=TRUE),
                   lower = mean-sd, 
                   upper = mean+sd
                   )
meangrowth

##    Community   N     mean       sd        se max    lower    upper
## 1 Individual  45 3.533333 1.778661 0.2651472   8 1.754672 5.311995
## 2      Group 112 3.357143 1.558800 0.1472928   7 1.798343 4.915943

meangrowth2 <- ddply(final, c("Fusion.Partners"), summarise,
                   N    = length(Polyps[!is.na(Polyps)]),
                   mean = mean(Polyps, na.rm=TRUE),
                   sd   = sd(Polyps, na.rm=TRUE),
                   se   = sd / sqrt(N), 
                   max = max(Polyps, na.rm=TRUE),
                   lower = mean-sd, 
                   upper = mean+sd
                   )
meangrowth2

##   Fusion.Partners  N     mean       sd        se max    lower    upper
## 1               0 71 3.267606 1.715124 0.2035478   8 1.552482 4.982730
## 2               1 47 3.659574 1.605379 0.2341687   7 2.054195 5.264954
## 3               2 23 3.043478 1.460954 0.3046300   7 1.582524 4.504433
## 4               3 16 3.812500 1.376893 0.3442232   6 2.435607 5.189393

meangrowth3 <- ddply(final, c("Richness"), summarise,
                   N    = length(Polyps[!is.na(Polyps)]),
                   mean = mean(Polyps, na.rm=TRUE),
                   sd   = sd(Polyps, na.rm=TRUE),
                   se   = sd / sqrt(N), 
                   max = max(Polyps, na.rm=TRUE),
                   lower = mean-sd, 
                   upper = mean+sd
                   )
meangrowth3

##   Richness  N     mean       sd        se max    lower    upper
## 1        1 86 3.290698 1.540560 0.1661228   8 1.750138 4.831257
## 2        2 16 3.937500 1.611159 0.4027898   7 2.326341 5.548659
## 3        3 40 3.525000 1.663908 0.2630870   7 1.861092 5.188908
## 4        4 15 3.200000 1.971222 0.5089672   7 1.228778 5.171222

#generate summary table 
grow_polyps <- ddply(final, c("Fusion.Type"), summarise,
                   N    = length(Polyps[!is.na(Polyps)]),
                   mean = mean(Polyps, na.rm=TRUE),
                   sd   = sd(Polyps, na.rm=TRUE),
                   se   = sd / sqrt(N)) 
                
grow_polyps

##   Fusion.Type  N     mean       sd        se
## 1       Fused 86 3.523256 1.539494 0.1660079
## 2  Individual 45 3.533333 1.778661 0.2651472
## 3     Unfused 26 2.807692 1.523660 0.2988142

Juveniles generated approximately 3-4 new polyps over the course of the grow-out period.

Plotting

Plot the model output of the effect of fusion partners and genotypic richness on polyp expansion.

eff.polyps1 <- predictorEffect("Fusion.Partners", polypmod1, xlevels=4, rug=TRUE) #set x axis
polyp.effects1<-plot(eff.polyps1,
                   lines=list(multiline=TRUE, #color lines
                              col=c("lightblue","mediumblue", "blue4", "darkblue"), 
                              lty=1), 
                   confint=list(style="bands", alpha=0.1), #set conf int style
                   lwd=4,
                   type="response", #set response scale
                   ylab=expression(bold("Total Polyp Expansion")), 
                   legend.position="right",
                   xlab=expression(bold("Fusion Partners")), 
                   main="",
                   lattice=list(key.args=list(space="right",
                                              border=FALSE, 
                                              title=expression(bold("Genotypic \nRichness")),
                                              cex=1, 
                                              cex.title=1)));polyp.effects1

4. Grow Out Period Figures

Combine figures for Grow-Out Period using cow plot and export to figures folder.

grow_out_figB<-plot_grid(surv.effects2b, labels = c(""), ncol=1, nrow=1, rel_heights= c(1), rel_widths = c(1), label_size = 20, label_y=1, align="h");grow_out_figB

ggsave(filename="Figures/grow_out_figB.png", plot=grow_out_figB, dpi=500, width=8, height=8, units="in")
ggsave(filename="Figures/grow_out_figB.pdf", plot=grow_out_figB, dpi=500, width=8, height=8, units="in")

grow_out_figs<-plot_grid(mean_growout_plot, surv.effects2b, labels = c("A", "B"), ncol=2, nrow=1, rel_heights= c(1,0.8), rel_widths = c(0.6,1), label_size = 20, label_y=1, align="v");grow_out_figs

## Warning: Graphs cannot be vertically aligned unless the axis parameter is set.
## Placing graphs unaligned.

ggsave(filename="Figures/grow_out_figs.png", plot=grow_out_figs, dpi=500, width=12, height=8, units="in")
ggsave(filename="Figures/grow_out_figs.pdf", plot=grow_out_figs, dpi=500, width=12, height=8, units="in")

5. Physical Data: Temperature Conditions

Temperature Conditions

Temperature data in this experiment were recorded every 15 minutes by a Neptune Apex System. Neptune probes were calibrated to a digital thermometer weekly.

Load in temperature data from Apex files. Display mean daily maximum temperatures in each treatment.

juvtemps<-read.csv("Data/ApexTemps.csv", header=T, na.strings="NA") #load data
#convert Date and Time to date/time format
juvtemps$Date <- as.POSIXct(juvtemps$Date, format="%m/%d/%y %H:%M")
juvtemps$DateDay <- as.POSIXct(juvtemps$DateDay, format="%m/%d/%y")

#convert to long format
juvtemps<-juvtemps %>% gather(Tank, Temperature, Tank17:Tank24)
juvtemps$Tank<-as.factor(juvtemps$Tank)

dailymax<-juvtemps[which(juvtemps$TimeDay == "930"),] #subset for daily maximum temperature
dailymax<-juvtemps[which(juvtemps$Phase == "Stress"),] #subset for stress period

dailymax$Treatment<-ifelse(dailymax$Tank %in% c("Tank17"), "Ambient", #assign tank treatments
                            ifelse(dailymax$Tank %in% c("Tank18"), "Ambient", 
                            ifelse(dailymax$Tank %in% c("Tank19"), "High",
                            ifelse(dailymax$Tank %in% c("Tank20"), "High",
                            ifelse(dailymax$Tank %in% c("Tank21"), "Ambient",
                            ifelse(dailymax$Tank %in% c("Tank22"), "High",
                            ifelse(dailymax$Tank %in% c("Tank23"), "Ambient",
                            ifelse(dailymax$Tank %in% c("Tank24"), "High",NA))))))))

maxtemps <- ddply(dailymax, c("Treatment"), summarise,
                   N    = length(Temperature[!is.na(Temperature)]),
                   mean = mean(Temperature, na.rm=TRUE),
                   sd   = sd(Temperature, na.rm=TRUE),
                   se   = sd / sqrt(N), 
                   max = max(Temperature, na.rm=TRUE),
                   lower = mean-sd, 
                   upper = mean+sd
                   )
maxtemps

juvtemps$Period<-ifelse(juvtemps$Phase %in% c("Stress"), "Stress", #assign stress period names
                            ifelse(juvtemps$Phase %in% c("Stress2"), "Stress", 
                            ifelse(juvtemps$Phase %in% c("Acclimation"), "Stress",
                            ifelse(juvtemps$Phase %in% c("Hurricane"), "Hurricane",
                            ifelse(juvtemps$Phase %in% c("Ambient"), "Grow-Out",NA)))))

Plot a timeseries of temperature data with tank temperature in colors (red=high, blue=ambient).

#subset temperature measurements

juvtempsHIGH<-subset(juvtemps, Period==c("Stress"), c(Date, DateDay, TimeDay, Phase, Tank, Temperature, Period))
juvtempsHIGH$Treatment<-ifelse(juvtempsHIGH$Tank %in% c("Tank17"), "Ambient",
                            ifelse(juvtempsHIGH$Tank %in% c("Tank18"), "Ambient", 
                            ifelse(juvtempsHIGH$Tank %in% c("Tank19"), "High",
                            ifelse(juvtempsHIGH$Tank %in% c("Tank20"), "High",
                            ifelse(juvtempsHIGH$Tank %in% c("Tank21"), "Ambient",
                            ifelse(juvtempsHIGH$Tank %in% c("Tank22"), "High",
                            ifelse(juvtempsHIGH$Tank %in% c("Tank23"), "Ambient",
                            ifelse(juvtempsHIGH$Tank %in% c("Tank24"), "High",NA))))))))

juvtempsHURR<-subset(juvtemps, Period==c("Hurricane"), c(Date, DateDay, TimeDay, Phase, Tank, Temperature, Period))
juvtempsHURR$Treatment<-ifelse(juvtempsHURR$Tank %in% c("Tank17"), "Ambient",
                            ifelse(juvtempsHURR$Tank %in% c("Tank18"), "Ambient", 
                            ifelse(juvtempsHURR$Tank %in% c("Tank19"), "High",
                            ifelse(juvtempsHURR$Tank %in% c("Tank20"), "High",
                            ifelse(juvtempsHURR$Tank %in% c("Tank21"), "Ambient",
                            ifelse(juvtempsHURR$Tank %in% c("Tank22"), "High",
                            ifelse(juvtempsHURR$Tank %in% c("Tank23"), "Ambient",
                            ifelse(juvtempsHURR$Tank %in% c("Tank24"), "High",NA))))))))

juvtempsGROW<-subset(juvtemps, Period==c("Grow-Out"), c(Date, DateDay, TimeDay, Phase, Tank, Temperature, Period))
juvtempsGROW$Treatment<-ifelse(juvtempsGROW$Tank %in% c("Tank17"), "Ambient",
                            ifelse(juvtempsGROW$Tank %in% c("Tank18"), "Ambient", 
                            ifelse(juvtempsGROW$Tank %in% c("Tank19"), "Ambient",
                            ifelse(juvtempsGROW$Tank %in% c("Tank20"), "Ambient",
                            ifelse(juvtempsGROW$Tank %in% c("Tank21"), "Ambient",
                            ifelse(juvtempsGROW$Tank %in% c("Tank22"), "Ambient",
                            ifelse(juvtempsGROW$Tank %in% c("Tank23"), "Ambient",
                            ifelse(juvtempsGROW$Tank %in% c("Tank24"), "Ambient",NA))))))))

juvtemps<-rbind(juvtempsHIGH, juvtempsGROW)
juvtemps<-rbind(juvtemps, juvtempsHURR)

TimeSeries1j<-ggplot(juvtemps, aes(x = Date, y = Temperature)) + 
  geom_line(aes(color = Treatment), size = 1) +
  scale_color_manual(values=c("blue", "red")) +
  ylab("Temperature (°C)")+
  xlab("Date")+
  ylim(26,32)+
  geom_vline(xintercept = as.POSIXct(as.Date(c("2018-07-31 10:00:00"))), linetype="solid", 
                color = "black", size=1)+
  geom_vline(xintercept = as.POSIXct(as.Date(c("2018-08-16 00:00:00"))), linetype="dashed", 
                color = "black", size=1)+
  geom_vline(xintercept = as.POSIXct(as.Date(c("2018-09-20 16:00:00"))), linetype="dotted", 
                color = "black", size=1)+
  theme_classic()+
  scale_x_datetime(limits = as.POSIXct(as.Date(c("2018-07-30 09:00:00", "2018-09-20 21:00:00")))) +
  theme(legend.position="none")+
  theme(text = element_text(size = 18, color="black"))+
  theme(axis.text = element_text(size=16, color="black"));TimeSeries1j

Solid line indicates start of experiment; dashed line indicates start of stress period; dotted line indicates end of experiment. Note the reduction in temperature in the middle of the stress test due to a system shut down as a result of a hurricane.

Plot by time of day by summarizing the mean temp at each time of day for each tank. Display the mean temperatures of each treatment during the Stress test phase (excluding the pause in treatment due to the hurricane system shutdown).

#subset juv temps for the stress period
juvtempsSTRESS<-juvtemps
juvtempsSTRESS$Period<-ifelse(juvtempsSTRESS$Phase %in% c("Stress"), "Stress",
                            ifelse(juvtempsSTRESS$Phase %in% c("Stress2"), "Stress", 
                            ifelse(juvtempsSTRESS$Phase %in% c("Acclimation"), "Acclimation",
                            ifelse(juvtempsSTRESS$Phase %in% c("Hurricane"), "Hurricane",
                            ifelse(juvtempsSTRESS$Phase %in% c("Ambient"), "Grow-Out",NA)))))

juvtempsSTRESS<-juvtempsSTRESS[which(juvtempsSTRESS$Phase == "Stress"),]


QuarterlyJ <- ddply(juvtempsSTRESS, c("TimeDay", "Treatment"), summarise,
                   N    = length(Temperature[!is.na(Temperature)]),
                   mean = mean(Temperature, na.rm=TRUE),
                   sd   = sd(Temperature, na.rm=TRUE),
                   se   = sd / sqrt(N), 
                   max = max(Temperature, na.rm=TRUE),
                   lower = mean-sd, 
                   upper = mean+sd
                   )

MeansJtemp <- ddply(juvtemps, c("Period", "Treatment"), summarise,
                   N    = length(Temperature[!is.na(Temperature)]),
                   mean = mean(Temperature, na.rm=TRUE),
                   sd   = sd(Temperature, na.rm=TRUE),
                   se   = sd / sqrt(N)
                   )

MeansJtemp2 <- ddply(juvtemps, c("Treatment"), summarise,
                   N    = length(Temperature[!is.na(Temperature)]),
                   mean = mean(Temperature, na.rm=TRUE),
                   max = max(Temperature, na.rm=TRUE),
                   sd   = sd(Temperature, na.rm=TRUE),
                   se   = sd / sqrt(N)
                   )
MeansJtemp2

##   Treatment     N     mean    max        sd          se
## 1   Ambient 30133 27.33872 29.315 0.4721296 0.002719819
## 2      High 13401 29.38521 31.700 1.5396924 0.013300420

Generate final figure of temperature treatments averaged over a 24-hour period.

JuvTempsDay<-ggplot(QuarterlyJ, aes(x = TimeDay, y = mean, color=Treatment, fill=Treatment)) + 
  geom_line(color="black", size = 1) +
  #facet_wrap(~Period) +
  scale_x_continuous(breaks=c(0, 720, 1440), labels=c("0:00", "12:00", "24:00")) +
  scale_color_manual(values = c("blue", "red")) +
  scale_fill_manual(values = c("blue", "red")) +
  theme_classic() +
  geom_ribbon(aes(ymin=QuarterlyJ$lower, ymax=QuarterlyJ$upper), linetype=2, alpha=0.4, color=NA) +
  ylab("Temperature (°C)") +
  xlab("Time of Day") +
  ylim(26,32)+
  theme(text = element_text(size = 18, color="black"))+
  theme(panel.margin = unit(2, "lines"))+
  theme(axis.text = element_text(size=16, color="black")); JuvTempsDay

Combine figures for temperature.

temp_figs<-plot_grid(TimeSeries1j, JuvTempsDay, labels = c("A", "B"), ncol=2, nrow=1, rel_heights= c(1,1), rel_widths = c(1,0.8), label_size = 20, label_y=1, align="h")

ggsave(filename="Figures/temp_figs.png", plot=temp_figs, dpi=500, width=14, height=6, units="in")

Temperature profiles during the stress period.

6. Stress Period: Survivorship

Load data

Survivorship is measured as binary response (alive=1, dead=0) on each individual spat in the experiment at multiple timepoints over the course of the stress period.

Due to mortality in juveniles in groups without fusion partners (see section 1) we removed unfused juveniles in groups from the analysis as there was low sample size of these groups. During the stress period, due to reductions in sample size with mortality, we no longer classify fusion by the number of fusion partners, but rather a binary classification of “fused” or “individual”. Similarly, we classified genotypic richness either as “single” genotypes or “multiple” genotypes.

Load data.

#load data
stress<-read.csv("Data/Stress_Responses_PostGenotyping.csv", header=TRUE, sep=",", na.strings="NA") #load data
stress$Diversity.Type<-relevel(stress$Diversity.Type, ref="Individual") #format columns
stress$Community<-relevel(stress$Community, ref="Individual")
stress$Polyps<-as.numeric(stress$Polyps)
stress$Fusion.Type<-relevel(stress$Fusion.Type, ref="Not Fused")
stress$Richness.Type<-relevel(stress$Richness.Type, ref="Single")

#remove non-fused juveniles in groups
stressdrop<-droplevels(stress[!stress$Fusion.Type == 'Not Fused',])

Subset separate timepoints for additional analyses.

final_stress<-stress[which(stress$Timepoint == "End"),]
final_stress_drop<-stressdrop[which(stressdrop$Timepoint == "End"),]
mid_stress_drop<-stressdrop[which(stressdrop$Timepoint == "S2D8"),]
group_stress<-stress[which(stress$Community == "Group"),]
group_final_stress<-group_stress[which(group_stress$Timepoint == "End"),]
group_d26_stress<-group_stress[which(group_stress$Timepoint == "S2D11"),]
d26_stress<-stress[which(stress$Timepoint == "S2D11"),]
d23_stress<-stress[which(stress$Timepoint == "S2D8"),]

Analyze survivorship by day, with a logistic binomial regression model.

Main effects: Day, temperature treatment, fusion, genotypic richness
Random effects: Parent colony, tank, slide

stressmod1<-glmer(Survivorship ~ Day + Treatment + Fusion.Type + Richness.Type + Day:Treatment + Day:Fusion.Type + Day:Richness.Type + Day:Treatment:Fusion.Type + Day:Treatment:Richness.Type + Treatment:Fusion.Type + Treatment:Richness.Type + (1|Parent.Site/Colony) + (1|Treatment:Tank) + (1|Slide), data=stressdrop, family=binomial) 

summary(stressmod1)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: Survivorship ~ Day + Treatment + Fusion.Type + Richness.Type +  
##     Day:Treatment + Day:Fusion.Type + Day:Richness.Type + Day:Treatment:Fusion.Type +  
##     Day:Treatment:Richness.Type + Treatment:Fusion.Type + Treatment:Richness.Type +  
##     (1 | Parent.Site/Colony) + (1 | Treatment:Tank) + (1 | Slide)
##    Data: stressdrop
## 
##      AIC      BIC   logLik deviance df.resid 
##    329.1    406.9   -148.5    297.1      939 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.5294  0.0000  0.0080  0.0553  3.6877 
## 
## Random effects:
##  Groups             Name        Variance  Std.Dev.
##  Slide              (Intercept) 7.056e+00 2.65637 
##  Colony:Parent.Site (Intercept) 1.966e+00 1.40198 
##  Treatment:Tank     (Intercept) 2.375e-03 0.04873 
##  Parent.Site        (Intercept) 3.306e-05 0.00575 
## Number of obs: 955, groups:  
## Slide, 81; Colony:Parent.Site, 20; Treatment:Tank, 9; Parent.Site, 6
## 
## Fixed effects:
##                                          Estimate Std. Error z value Pr(>|z|)
## (Intercept)                              14.66027    6.10385   2.402  0.01631
## Day                                      -0.39207    0.19329  -2.028  0.04251
## TreatmentHigh                            56.33441   19.78098   2.848  0.00440
## Fusion.TypeIndividual                    -5.31708    6.35453  -0.837  0.40274
## Richness.TypeMultiple                     0.11981    7.18623   0.017  0.98670
## Day:TreatmentHigh                        -2.49527    0.81020  -3.080  0.00207
## Day:Fusion.TypeIndividual                 0.27467    0.20584   1.334  0.18207
## Day:Richness.TypeMultiple                 0.02826    0.23089   0.122  0.90259
## TreatmentHigh:Fusion.TypeIndividual     -55.35582   19.90805  -2.781  0.00543
## TreatmentHigh:Richness.TypeMultiple     -53.24767   20.43721  -2.605  0.00918
## Day:TreatmentHigh:Fusion.TypeIndividual   2.11216    0.78811   2.680  0.00736
## Day:TreatmentHigh:Richness.TypeMultiple   2.17115    0.81460   2.665  0.00769
##                                           
## (Intercept)                             * 
## Day                                     * 
## TreatmentHigh                           **
## Fusion.TypeIndividual                     
## Richness.TypeMultiple                     
## Day:TreatmentHigh                       **
## Day:Fusion.TypeIndividual                 
## Day:Richness.TypeMultiple                 
## TreatmentHigh:Fusion.TypeIndividual     **
## TreatmentHigh:Richness.TypeMultiple     **
## Day:TreatmentHigh:Fusion.TypeIndividual **
## Day:TreatmentHigh:Richness.TypeMultiple **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) Day    TrtmnH Fsn.TI Rch.TM Dy:TrH D:F.TI D:R.TM TH:F.T
## Day         -0.956                                                        
## TreatmntHgh -0.210  0.250                                                 
## Fsn.TypIndv -0.917  0.899  0.322                                          
## Rchnss.TypM -0.814  0.796  0.285  0.796                                   
## Dy:TrtmntHg  0.131 -0.194 -0.993 -0.253 -0.223                            
## Dy:Fsn.TypI  0.879 -0.930 -0.286 -0.947 -0.752  0.237                     
## Dy:Rchns.TM  0.783 -0.829 -0.269 -0.763 -0.961  0.225  0.782              
## TrtmnH:F.TI  0.209 -0.249 -0.992 -0.351 -0.282  0.984  0.314  0.267       
## TrtmnH:R.TM  0.201 -0.241 -0.971 -0.311 -0.375  0.965  0.278  0.358  0.963
## Dy:TrH:F.TI -0.146  0.205  0.984  0.279  0.224 -0.990 -0.273 -0.227 -0.991
## Dy:TrH:R.TM -0.134  0.195  0.971  0.248  0.296 -0.978 -0.235 -0.304 -0.963
##             TH:R.T D:TH:F
## Day                      
## TreatmntHgh              
## Fsn.TypIndv              
## Rchnss.TypM              
## Dy:TrtmntHg              
## Dy:Fsn.TypI              
## Dy:Rchns.TM              
## TrtmnH:F.TI              
## TrtmnH:R.TM              
## Dy:TrH:F.TI -0.956       
## Dy:TrH:R.TM -0.991  0.970
## convergence code: 0
## Model failed to converge with max|grad| = 0.380499 (tol = 0.001, component 1)
## Model is nearly unidentifiable: large eigenvalue ratio
##  - Rescale variables?

Anova(stressmod1, type="II")

## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: Survivorship
##                               Chisq Df Pr(>Chisq)    
## Day                         22.2735  1  2.365e-06 ***
## Treatment                   11.8732  1  0.0005695 ***
## Fusion.Type                  0.3816  1  0.5367478    
## Richness.Type                0.2412  1  0.6233337    
## Day:Treatment               11.4821  1  0.0007027 ***
## Day:Fusion.Type              4.6173  1  0.0316509 *  
## Day:Richness.Type            0.9590  1  0.3274329    
## Treatment:Fusion.Type        0.8593  1  0.3539263    
## Treatment:Richness.Type      0.0686  1  0.7933137    
## Day:Treatment:Fusion.Type    7.1827  1  0.0073612 ** 
## Day:Treatment:Richness.Type  7.1038  1  0.0076922 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

There are significant three-way interactions between day, temperature, and fusion as well as day, temperature, and genotypic richness.

Check for overdisperson. No evidence of overdispersion.

dispersion_glmer(stressmod1) #no overdispersion

Plot with fusion and richness separately

First plot the effect of fusion on survival.

#generate model for effect plot  
stressdrop$group<-paste(stressdrop$Treatment, "-", stressdrop$Fusion.Type) #set group names as combinations of temperature and fusion to plot on one plot

stressmod2<-glmer(Survivorship ~ Day * group + (1|Parent.Site/Colony) + (1|Treatment:Tank) + (1|Slide), data=stressdrop, family=binomial, link=logit) #run with group names

eff.stress2 <- predictorEffect(c("Day"), stressmod2) #set x axis

stress.effectsFUSE<-plot(eff.stress2,
                   lines=list(multiline=TRUE, #color lines
                              col=c("blue", "blue", "red", "red"), 
                              lty=c(2,1,2,1)), 
                   confint=list(style="bands", alpha=0.1), #set conf int
                   lwd=4,
                   axes=list(y=list(lim=c(0, 1.1))),
                   type="response", #set response scale
                   ylab=expression(bold("Prob(Survivorship)")), 
                   legend.position="top",
                   xlab=expression(bold("Day")), 
                   main="",
                   lattice=list(key.args=list(space="right",
                                              border=FALSE, 
                                              title=expression(bold("Temperature - Fusion")),
                                              cex=1, 
                                              cex.title=1)));stress.effectsFUSE

Obtain effect summary for effects of fusion on survival.

summary(effect(c("Day", "group"), xlevels=50, stressmod2))

The curves for high temperature meet 50% at the following timepoints:

High individual = ~ 20.4 days
High fused = ~ 25.5 days

50% mortality is approximately 5 days later for fused corals than for individuals in this stress test.

Next, plot the effect of genotypic richness within the fused corals.

stressdrop$group<-paste(stressdrop$Treatment, "-", stressdrop$Richness.Type) #set group name for single panel plot

stressmod3<-glmer(Survivorship ~ Day * group + (1|Parent.Site/Colony) + (1|Treatment:Tank) + (1|Slide), data=stressdrop, family=binomial, link=logit, subset=c(Fusion.Type=="Fused")) #run model with group

eff.stress3 <- predictorEffect(c("Day"), stressmod3) #set x axis
stress.effectsRICHNESS<-plot(eff.stress3,
                   lines=list(multiline=TRUE, #set line color
                              col=c("blue", "blue", "red", "red"), 
                              lty=c(2,1,2,1)), 
                   confint=list(style="bands", alpha=0.1), #set conf int
                   lwd=4,
                   axes=list(y=list(lim=c(0, 1.1))),
                   type="response", #set resposne scale
                   ylab=expression(bold("Prob(Survivorship)")), 
                   legend.position="top",
                   xlab=expression(bold("Day")), 
                   main="",
                   lattice=list(key.args=list(space="right",
                                              border=FALSE, 
                                              title=expression(bold("Temperature - Richness")),
                                              cex=1, 
                                              cex.title=1)));stress.effectsRICHNESS

Obtain effect summary for effects of richness.

summary(effect(c("Day", "group"), xlevels=50, stressmod3))

The curves for high temperature meet 50% at the following timepoints:

High multiple = ~ 26 days
High single = ~ 24.5 days

50% mortality (LT50) is approximately 1.5 days later for multiple genotype fusions than for single genotype fusions in this stress test.

7. Stress Period: Fusion Persistence

During the stress period, we monitored the persistence of fusion as a binary response (fused = 1, not fused = 0).

Analysis

Analyze fusion persitence as an effect of treatment and genotypic richness with a logistic binomial regression model at the end of the stress period.

Main effects: Temperature treatment, genotypic richness
Random effects: Parent colony nested within parent collection site, tank, slide

Display the model summary and anova tables for use of p-values.

stressfuse1<-glmer(Fusion ~ Treatment * Richness.Type + (1|Parent.Site/Colony) + (1|Treatment:Tank) + (1|Slide), data=final_stress_drop, family=binomial, subset=c(Fusion.Type=="Fused")) 

summary(stressfuse1)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: Fusion ~ Treatment * Richness.Type + (1 | Parent.Site/Colony) +  
##     (1 | Treatment:Tank) + (1 | Slide)
##    Data: final_stress_drop
##  Subset: c(Fusion.Type == "Fused")
## 
##      AIC      BIC   logLik deviance df.resid 
##     95.2    115.3    -39.6     79.2       84 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -1.79687  0.05947  0.07494  0.15442  0.89648 
## 
## Random effects:
##  Groups             Name        Variance  Std.Dev. 
##  Slide              (Intercept) 1.741e+01 4.173e+00
##  Colony:Parent.Site (Intercept) 2.833e-09 5.322e-05
##  Treatment:Tank     (Intercept) 3.526e-08 1.878e-04
##  Parent.Site        (Intercept) 0.000e+00 0.000e+00
## Number of obs: 92, groups:  
## Slide, 36; Colony:Parent.Site, 20; Treatment:Tank, 8; Parent.Site, 6
## 
## Fixed effects:
##                                     Estimate Std. Error z value Pr(>|z|)  
## (Intercept)                           5.4603     2.9581   1.846   0.0649 .
## TreatmentHigh                        -2.2330     2.8363  -0.787   0.4311  
## Richness.TypeMultiple                -0.4727     2.5539  -0.185   0.8532  
## TreatmentHigh:Richness.TypeMultiple  -0.6406     3.6936  -0.173   0.8623  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) TrtmnH Rch.TM
## TreatmntHgh -0.488              
## Rchnss.TypM -0.552  0.629       
## TrtmnH:R.TM  0.191 -0.787 -0.710
## convergence code: 0
## boundary (singular) fit: see ?isSingular

Anova(stressfuse1, type="II") 

## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: Fusion
##                          Chisq Df Pr(>Chisq)
## Treatment               2.2446  1     0.1341
## Richness.Type           0.1914  1     0.6617
## Treatment:Richness.Type 0.0301  1     0.8623

There is no effect of treatment or genotypic richness on fusion persistence.

Check for overdispersion.

dispersion_glmer(stressfuse1) #no overdispersion

## [1] 0.6815656

Display proportion of juveniles fused by the end of the stress period. Greater than 70% of juveniles remained fused during this period.

meanstressfuse <- ddply(final_stress_drop, c("Treatment", "Richness.Type"), summarise,
                   N    = length(Fusion[!is.na(Fusion)]),
                   mean = mean(Fusion, na.rm=TRUE),
                   sd   = sd(Fusion, na.rm=TRUE),
                   se   = sd / sqrt(N), 
                   max = max(Fusion, na.rm=TRUE),
                   lower = mean-sd, 
                   upper = mean+sd
                   )
meanstressfuse

##   Treatment Richness.Type  N      mean        sd         se max     lower
## 1   Ambient        Single 12 0.9166667 0.2886751 0.08333333   1 0.6279915
## 2   Ambient      Multiple 29 0.8620690 0.3509312 0.06516629   1 0.5111378
## 3      High        Single 19 0.7368421 0.4524139 0.10379087   1 0.2844282
## 4      High      Multiple 32 0.7187500 0.4568034 0.08075220   1 0.2619466
##      upper
## 1 1.205342
## 2 1.213000
## 3 1.189256
## 4 1.175553

8. Stress Period: Polyp Expansion (Growth)

Analysis

Analyze juvenile growth by polyp expansionwith a poisson mixed effects model.

Analyze growth as polyp expansion at day 23 of 32 of the stress test. Due to high mortality, this is the last time point with suffficient sample size to measure growth. Growth is measured as the number of polyps per juvenile. As there was no difference in groups at the end of the grow-out period, we measure growth during the stress period as total number of polyps per juvenile.

Display mean number of polyps per juvenile at day 23.

meanstressgrowth <- ddply(mid_stress_drop, c("Treatment", "Fusion.Type", "Richness.Type"), summarise,
                   N    = length(Polyps[!is.na(Polyps)]),
                   mean = mean(Polyps, na.rm=TRUE),
                   sd   = sd(Polyps, na.rm=TRUE),
                   se   = sd / sqrt(N), 
                   max = max(Polyps, na.rm=TRUE),
                   lower = mean-sd, 
                   upper = mean+sd
                   )
meanstressgrowth

##   Treatment Fusion.Type Richness.Type  N     mean       sd        se max
## 1   Ambient       Fused        Single 12 6.416667 1.564279 0.4515685   8
## 2   Ambient       Fused      Multiple 29 6.827586 1.909917 0.3546627  10
## 3   Ambient  Individual        Single 20 7.000000 2.294157 0.5129892  11
## 4      High       Fused        Single 16 5.375000 1.668333 0.4170831   8
## 5      High       Fused      Multiple 29 5.965517 1.917639 0.3560967  10
## 6      High  Individual        Single 16 6.250000 1.570563 0.3926406   9
##      lower    upper
## 1 4.852387 7.980946
## 2 4.917669 8.737503
## 3 4.705843 9.294157
## 4 3.706667 7.043333
## 5 4.047878 7.883156
## 6 4.679437 7.820563

Analyze the effect of temperature, fusion, and genotypic richness on polyp expansion at day 23. Analyze with a mixed effect poisson model.

Main effects: Temperature treatment, genotypic richness, fusion
Random effects: Parent colony nested within parent collection site, tank, slide

Display model summary and anova table with p-values.

stresspolypmod1<-glmer(Polyps ~ Treatment + Richness.Type + Fusion.Type + Treatment:Richness.Type + Treatment:Fusion.Type + (1|Parent.Site/Colony) + (1|Treatment:Tank) + (1|Slide), data=stressdrop, family = poisson, subset=c(Timepoint=="S2D8"))

summary(stresspolypmod1)

## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: poisson  ( log )
## Formula: 
## Polyps ~ Treatment + Richness.Type + Fusion.Type + Treatment:Richness.Type +  
##     Treatment:Fusion.Type + (1 | Parent.Site/Colony) + (1 | Treatment:Tank) +  
##     (1 | Slide)
##    Data: stressdrop
##  Subset: c(Timepoint == "S2D8")
## 
##      AIC      BIC   logLik deviance df.resid 
##    537.1    565.2   -258.6    517.1      112 
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.12141 -0.41736  0.00521  0.46001  1.49721 
## 
## Random effects:
##  Groups             Name        Variance Std.Dev.
##  Slide              (Intercept) 0.00000  0.00000 
##  Colony:Parent.Site (Intercept) 0.00604  0.07771 
##  Treatment:Tank     (Intercept) 0.00000  0.00000 
##  Parent.Site        (Intercept) 0.00000  0.00000 
## Number of obs: 122, groups:  
## Slide, 72; Colony:Parent.Site, 20; Treatment:Tank, 8; Parent.Site, 6
## 
## Fixed effects:
##                                     Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                          1.88987    0.12770  14.799   <2e-16 ***
## TreatmentHigh                       -0.20863    0.17044  -1.224    0.221    
## Richness.TypeMultiple                0.02786    0.14665   0.190    0.849    
## Fusion.TypeIndividual                0.04796    0.15485   0.310    0.757    
## TreatmentHigh:Richness.TypeMultiple  0.06667    0.19937   0.334    0.738    
## TreatmentHigh:Fusion.TypeIndividual  0.10252    0.21773   0.471    0.638    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) TrtmnH Rch.TM Fsn.TI TH:R.T
## TreatmntHgh -0.750                            
## Rchnss.TypM -0.861  0.658                     
## Fsn.TypIndv -0.818  0.628  0.725              
## TrtmnH:R.TM  0.631 -0.850 -0.728 -0.525       
## TrtmnH:F.TI  0.600 -0.793 -0.529 -0.721  0.672
## convergence code: 0
## boundary (singular) fit: see ?isSingular

Anova(stresspolypmod1, type="II")

## Analysis of Deviance Table (Type II Wald chisquare tests)
## 
## Response: Polyps
##                          Chisq Df Pr(>Chisq)  
## Treatment               3.7610  1    0.05246 .
## Richness.Type           0.3997  1    0.52724  
## Fusion.Type             0.8766  1    0.34913  
## Treatment:Richness.Type 0.1118  1    0.73807  
## Treatment:Fusion.Type   0.2217  1    0.63775  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

There is no effect of treatment, fusion or genotypic richness on growth. Treatment is near significant at p=0.0525.

Assess normality assumption. Meets assumption.

qqPlot(residuals(stresspolypmod1))

## 752 750 
## 103 101

Plotting

Generate a plot of polyp expansion by temperature treatment to view effect.

stresspolypmod2<-glmer(Polyps ~ Treatment + (1|Parent.Site/Colony) + (1|Treatment:Tank) + (1|Slide), data=stressdrop, family = poisson, subset=c(Timepoint=="S2D8"))

eff.stress.polyps1 <- predictorEffect("Treatment", stresspolypmod2)
polyps.stress.effects1<-plot(eff.stress.polyps1,
                   lines=list(multiline=TRUE, 
                              col=c("black"), 
                              lty=1), 
                   confint=list(style="bars", width=4, col="black"), 
                   lwd=4,
                   axes=list(y=list(lim=c(5,8))),
                   type="response", 
                   ylab=expression(bold("Polyp Expansion")), 
                   legend.position="right",
                   xlab=expression(bold("Temperature")), 
                   main=""); polyps.stress.effects1

9. Stress Period Figures

Combine figures for Stress Period:

stress_figs<-plot_grid(stress.effectsFUSE, stress.effectsRICHNESS, labels = c("A", "B"), ncol=1, nrow=2, rel_heights= c(1,1), rel_widths = c(1,1), label_size = 20, label_y=1, align="h")

ggsave(filename="Figures/stress_figs.png", plot=stress_figs, dpi=500, width=12, height=12, units="in")

ggsave(filename="Figures/stress_figs.pdf", plot=stress_figs, dpi=500, width=12, height=12, units="in")

Effects of fusion and genotypic richness on juvenile coral survival during the stress period.